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Abstract 

We study the dynamics of the critical two-dimensional fully-frustrated 
Ising model by means of Monte Carlo simulations. The dynamical exponent is 
estimated at equilibrium and is shown to be compatible with the value — 2. 
In a second step, the system is prepared in the paramagnetic phase and then 
quenched at its critical temperature Tc = 0. Numerical evidences for the ex- 
istence of logarithmic corrections in the aging regime are presented. These 
corrections may be related to the topological defects observed in other fully- 
frustrated models. The autocorrelation exponent is estimated to be A = d as 
for the Ising chain quenched at Tc = 0. 

Introduction 

Out-of-equilibrium statistical physics is a very active field of research. Among the 
most studied topics, the slow evolution of glasses remains a challenging problem. 
During a quench in the glass phase, the system is trapped in an infinite succession 
of metastable states and never reaches equilibrium. Time-translation invariance 
is broken and the fiuctuation-dissipation theorem does not hold anymore [1, 2, 3]. 
Frustration and randomness are two key ingredients in the behavior of glasses. In- 
terestingly, neither frustration nor randomness are necessary conditions for aging. 
Indeed it is known for a decade that homogeneous ferromagnets can also display 
aging. When they are prepared in their paramagnetic phase and then quenched be- 
low their critical temperature Tc for example, ferromagnets age due to metastable 
states consisting in ferromagnetic domains in competition. The domain walls be- 
tween them cannot be eliminated in a finite time [4]. The divergence of the re- 
laxation time provokes the suppression of the exponential relaxation and two-time 
observables, for instance the autocorrelation function C(t, s), decays as a power-law 
of the scaling variable t/s. A more detailed presentation will be given in section 3. 

In this work, we are interested in the intermediate situation of frustrated sys- 
tems but without disorder. The model considered is the critical two-dimensional 
fully- frustrated Ising model (FFIM) . The equilibrium properties of the FFIM will be 
presented in section 1. The FFIM belongs to the same universality class as the anti- 
ferromagnetic Ising model on a triangular lattice (AFIT). As a consequence, one 
may assume that these two models display the same behavior out-of-equilibrium. 
This is indeed claimed in a note in Ref. [5]. However, controversial results are 
found in the literature about the AFIT. On the one hand, numerical evidences sug- 
gesting the existence of topological defects in anti-ferromagnetic Ising and Potts 
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Figure 1: Coupling configurations discussed in the text: Zig-Zag (left) and Piled-up 
Domino (right) configurations. Thin lines correspond to ferromagnetic couplings 
and thick ones to anti-ferromagnetic couplings. 



models have been given [6, 7]. Out-of-equilibrium, these topological defects behave 
as impurities and pin the domain walls. The motion of the latter is thus slowed 
down by a logarithmic factor which is then observable in dynamical quantities [4]. 
Such a mechanism is also encountered in pure systems where topological defects 
are present too. The paradigmatic model in this case is the two-dimensional XY 
model for which it has been shown that vortices manifest themselves by logarith- 
mic corrections in the long-time behavior of the dynamical observables [4, 8] . Such 
logarithmic corrections have been observed [7] in the AFIT. Like in the XY model, 
the dynamical exponent is predicted to be Zc = 2 in the AFIT since it exists an 
exact mapping onto a Gaussian Solid-on-Solid model [9]. On the other hand, a re- 
cent Monte Carlo simulation of the anti-ferromagnetic Ising model failed to observe 
such logarithmic corrections [5]. In contradistinction to the previous scenario, the 
dynamical exponent was estimated to be Zc — 2.33. The spin-spin autocorrelation 
was observed to decay as a power-law with an exponent Xc/zc — 0.86. 

Our goal is to investigate the existence of logarithmic corrections in the FFIM. 
To that purpose, we will estimate the dynamical exponent Zc at equilibrium in 
section 2. Remind that the two scenarii give quite different values. In section 3, we 
will study the aging properties and search for the signature of topological defects 
when the system is quenched from the paramagnetic phase. The ratio X/z will be 
estimated and compared to the different scenarii. 



1 The Fully-Frustrated Ising Model (FFIM) 

We consider the classical two-dimensional fully-frustrated Ising model (FFIM) de- 
fined by the Hamiltonian 

= - X] + i-iy'''^'^^crx,yax^y+i , a^^y G {-1, 1}. (1) 

The function f{x, y) is chosen such that each plaquette has an odd number of anti- 
ferromagnetic couplings. As a consequence there exists no spin configuration for 
which all bonds can be satisfied. The most studied coupling configurations corre- 
spond to f{x,y) = X + y (Zig-Zag) and f{x,y) = x (Piled-up Domino). They are 
depicted on Figure 1. At equilibrium, these two coupling configurations arc equiva- 
lent since a transformation leaving the partition function unchanged maps one onto 
the other [10]. In this study, we will consider the Zig-Zag coupling configuration. 



Spin configurations in the ground state can be obtained by minimization of the 
energy of each plaquette independently. Since by construction of the model it is 
not possible to satisiy all the bonds in a plaquette, the possible spin configurations 
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Figure 2: Possible spin configurations on a plaquette in the ground state. The 
dashed hues correspond to the domain walls between the different ferromagnetic 
phases. 
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Figure 3: Examples of spin configurations in the ground state of the FFIM for the 
Zig-Zag (left) and the Piled- up Domino coupling configuration (right). 



in the ground state are those with only one unsatisfied bond. They are shown on 
Figure 2. The ground-state is then built by juxtaposition of these plaquettes. Even 
when three out of the four spins of the plaquette are fixed by neighboring plaquettes, 
there may exist several possibilities to choose the fourth spin (the two last plaquettes 
of Figure 2 give an example of such a situation) . As a consequence, the ground state 
is highly degenerate and the entropy per site is finite at T = [11]. One can describe 
the spin configurations at T = as made of ferromagnetic domains separated by 
domain walls. The possible plaquettes in the ground state forbid that more than 
one domain wall goes through a plaquette which means that domain walls cannot 
intersect each other. Moreover, when going through a plaquette, the domain wall 
must cross the anti-ferromagnetic bond and one of the three ferromagnetic bonds. 
Figure 3 shows examples of spin configurations in the ground state. In the case of 
the Zig-Zag coupling configuration, domain walls cannot form overhangs or loops. 
They must start and end at two different boundaries of the system. For the Piled- up 
Domino coupling configuration, domain walls can form overhangs and loops. They 
can thus enclose finite clusters that can be as small as one single spin. 



The equilibrium properties of the FFIM have been determined by analytical di- 
agonalization of the transfer matrix [12]. The ferromagnetic order is destroyed at 
any temperature because the domain walls discussed above have no energy. Nev- 
ertheless, the system is critical at Tc = with spin-spin correlation functions de- 
caying algebraically with an exponent = 1/2 [13]. As already mentioned in the 
introduction, the FFIM belongs to the same universality class as the homogeneous 
anti- ferromagnetic Ising model on a triangular lattice [14] (AFIT). 
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2 Dynamical exponent of the FFIM 



We now consider the dynamical properties of the FFIM at equilibrium at its critical 

temperature = 0. The system is initially thermalized and then evolves with the 
Glauber dynamics [15]. For an inverse temperature j3 = l/ksT, the probability for 
a spin to be flipped is .^^^.^ae where AE is the variation of energy when the spin is 
flipped. In our case, the temperature is T = which means that a spin flip is always 
accepted if the energy is lowered, accepted with a probability 1/2 if the energy does 
not change, and always rejected otherwise. Moreover, we start from an equilibrium 
state at r = 0; i.e. with the minimal value of the energy. As a consequence, the 
algorithm becomes particularly simple: a spin is randomly chosen and is flipped if 
the total energy does not change ^. 

In contradistinction to homogeneous ferromagnets, the Glauber dynamics of the 
FFIM is not frozen at T = 0. In the case of the Zig-Zag coupling configuration, some 
single-spin flips, as for instance the formation of a dot on a flat domain wall, does 
not change the total enc;rgy and are thus allow(xi. In the; case; of the Pile-up Domino 
coupling configuration, one can also flip any spin which lies on an anti-ferromagnetic 
bond in the bulk of a ferromagnetic domain. This process creates a new single-spin 
domain as the one depicted on Figure 3. These single-spin domains can be flipped 
without any energy change. As a consequence, they behave as isolated Ising spins. 
This situation was already encountered in the case of the AFIT where these spins 
were called loose spins [13]. 

The dynamical exponent Zc of the FFIM at T = is not known. Since the FFIM 
belongs to the same universality class at equilibrium as the AFIT, one may assume 
that they share the same dynamical exponents Zc- As already mentioned in the 
introduction, different values can be found in the literature. According to Moore 
et al. [7] who claim the existence of topological defects, the dynamical exponent 
should be Zc = 2. On the other hand, the numerical estimate Zc — 2.33 has been 
determined in more recent Monte Carlo simulations [5]. 

To measure the critical exponent in the FFIM, we computed the spin-spin au- 
tocorrelation function at equilibrium at Tc = 0. The lattice size is 192 x 192 and 
we used periodic boundary conditions in both directions. The system is flrst ther- 
malized using the Kandel-Ben Av-Domany cluster algorithm [16]. We then use the 
Glauber dynamics and measure (Ti{s) and (Ti(t). The simulation is repeated 1000 
times. The mean spin-spin autocorrelation is finally estimated as 

Ce^Xt,s) = ^Y.{ai{t)ai{s)) (2) 

i 

where (. . .) denotes the average over the 1000 histories of the system. As shown on 
Figure 4, the correlation function decays as a power-law for sufficiently large sepa- 
ration time {t — s > 10). One indeed expects that the divergence of the relaxation 
time at the critical point leads to an algebraic decay of the autocorrelation function 
C(t, s). Moreover, time-translation invariance at equilibrium imposes a dependence 
on t — s. Finally, the order parameter having an anomalous dimension P/v, one 
expects 

Ce^.it,s) ~ (3) 
t-s>l 

The interpolation of the data (Figure 4) with this expression leads to the estimate 
2fi/uzc = 0.2495(11). Since r? = 2l3/v = 1/2 {d = 2), our estimate for the dy- 
namical exponent is Zc ~ 2.004(9). This value is compatible with Zc = 2, i.e. the 
one expected in the scenario involving topological defects [7] and contradicts the 

Changing the probabiUty of a spin flip from 1/2 to 1 leads here only to a factor of 2 in the 
autocorrelation time. 
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Figure 4: Equilibrium spin-spin autocorrelation function C{t,s) of the FFIM at 
Tc ~ with the Zig-Zag coupling configuration versus t — s. The six curves cor- 
responding to s = 10, 20, 40, 160, 640 and s = 1280 cannot be distinguished as 
expected at equilibrium. A very nice power-law is observed and allows for a precise 
determination of 2(3 /vZc- 

observation of a sub-diffusive growth with z ~ 2.33 observed by Kim et al. [5] with 
random initial configurations in the AFIT. Note that we also measured the dy- 
namical exponent Zc during a quench by considering the decay of autocorrelation 
functions in the quasi-equilibrium regime. We obtained again a value compatible 
with Zc = 2. 



We are now interested in the out-of-equilibrium dynamics of the FFIM. The system 
is initially prepared at T ^ -|-oo, i.e. spins are random. It it then quenched at the 
critical temperature Tc = 0. We use the Glauber dynamics at Tc — 0: a spin flip 
is always accepted when it leads to a decrease of the total energy {AE < 0) and 
only with a probability 1/2 if AE = 0. A typical spin configuration is presented 
on Figure 5. As expected, ferromagnetic domains grow during the quench and 
the domain walls tend to the stable configurations discussed in section 1. On the 
snapshot (Figure 5), most of the domains are already in the ground state. A few 
others do not span the entire systcnn and are thus unstable. One expects that thciy 
will disappear by shrinking or moving to one of the boundaries. In the following, we 
will show that the aging of the FFIM can be analyzed using the same assumptions 
as for homogeneous ferromagnets. 



In the case of homogeneous ferromagnets quenched below their critical tempera- 
ture, domains grow with a characteristic length L{t) behaving in time as L[t) ~ t^^^ 



3 Aging of the FFIM 
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Figure 5: Typical spin configuration of the FFIM with the Zig-Zag coupling config- 
uration when initially prepared at T ^ +00 and then quenched at Tc = 0. 

where z = 2 is the dynamical exponent [17]. The two-time autocorrelation function 
of the local order parameter can be decomposed as 



C{t,s) = C,t.{t-s) + C,^.{t/s) 



(4) 



where Cst. (t — s) is due to reversible processes occurring inside the domains while 
Cag. comes from the irreversible motion and annihilation of domain walls. This 
decomposition is clearly seen for example in the case of the spherical model [18]. 
Usually, the first term falls down rapidly and can be neglected. The self-similarity 
of the domain growth implies that two-time observables, say C(f, .s), depend only 
on L{t)/L{s) and thus on t/s when 1 ^ s ~ i.e. in the aging regime [4]. The 
autocorrelation function C{t, s) behaves as [19, 20]: 



(5) 



with the scaling variable t/s. In the limit t ^ s, the scaling function fc{x) decays 
algebraically as x~^l^ . This defines a new exponent, the autocorrelation exponent 
A [21]. At the critical point, similar conclusions can be drawn from the assumption 
that the correlation length grows during a quench as ^ ^ t^^^'= . The autocorrelation 
function C{t,s) displays a scaling behavior analogous to (5): 



C'ag. (i) s) 



"gc 



(6) 



where gci^) is expected to decay asymptotically as x ^"1^". Dynamical and au- 
tocorrelation exponents and Ac take values different from below T,,. Because 
C(s,s) = (M2(s)) ~ ^-2/3/"' ^ g-'^HIvz,^ ti^g exponent is equal to 2filvz^. The 
expression (6) works for a broad range of models for example the spherical model [18] 
or the Ising chain [22] that can be both treated analytically, the 2D Ising model 
by Monte Carlo simulations [18] or the 0(n) model analyzed by Renormalization 
Group technics [23]. In the case of the XY model, the motion of the domains is 
slowed down by the vortices that behave as impurities pinning the domain walls. 
The correlation length grows as (tj lnt)^/^= so that the autocorrelation function is 
modified by the replacement of the scaling variable tjs by f In s/s In t [24, 4]: 



Cag.(t,s) ~ s "'^'gc 



tins 
slnt 



(7) 
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We will analyze the aging of the FFIM in two steps. First, we will determine 
whether topological defects exist in the FFIM and affect the dynamics as claimed by 
some authors in the case of the AFIT. To that purpose, we will determine which one 
of the two scenarii (6) or (7) is the more compatible with our data, i.e. which scaling 
variable, t/s or tins/ slnt, leads to a collapse of the scaling function s°"=C{t,s) for 
different waiting times s. Then, we will estimate the autocorrelation exponent A 
using the appropriate scaling variable. For these simulations, the lattice size is again 
L = 192 with periodic boundary conditions but the data have been averaged over 
50,000 different histories of the system. 

3.1 Determination of the scaling vEiriable 

Our numerical estimation of the scaling fimction .s°'=C(f, s) is plotted with respect 
to t/s on Figure 6 and to tins/slut on Figure 7. We used the exponent = 
0.2495(11) obtained in section 2. No collapse of the scaling function s°"'C{t,s) is 
observed with the scaling variable t/s (Figure 6). The last two curves (s = 640 
and 1280) are however closer to each other. Our waiting times s are too small to 
check if this is the sign of an accumulation of the curves for much larger waiting 
times. In contradistinction, a very good collapse for large waiting times s > 160 can 
be seen in the inset of Figure 7 when the data are plotted versus t Ins/ sin t. The 
difference between the different curves is of the order of 2.10"^ where the correlation 
itself is of order of 0.2. The difference between the different curves is of the order 
of magnitude of the statistical fluctuations of the data. For these reasons, it is 
unlikely that this collapse for the three waiting times s = 160, s = 640 and 1280 be 
accidental and disappear for larger waiting times. We have also tried more complex 
scaling variables as a; = t{\ns)'^ /s{\nt)^ or x = Hn^/sln^ but they do not lead 
to a better collapse. 

Note that the observation of these logarithmic corrections has been made pos- 
sible by the use of the two times t and s. It is well known that it is very dif- 
ficult to distinguish these logarithmic corrections if ones studies only the asymp- 
totic dependence of C{t,s) on t, i.e. C{t,s) ^ {t/ \nt)~-^/^ . Moreover, if such 
corrections are observed, they may also result from corrections to scaling, i.e. 
C{t,s) ~ t--^/^(l + at-^ + ...) with w small. The relations (6) or (7) imposes 
a much stronger constraint than the asymptotic law C{t.s) ^ {t/ \Tit)~^/^ . We are 
looking for logarithmic factors in the scaling variable. The analysis is thus not sen- 
sitive to scaling corrections to the dominant algebraic decay. These considerations 
may explain why these logarithmic corrections, although observed by [7] for the 
AFIT have not be seen in more recent Monte Carlo simulations [5]. 



3.2 Determination of the exponent 

The previous analysis relies strongly on the scaling form (5). Our conclusion may 
be criticized by arguing that corrections to the prefactor s~"= should have been 
taken into account. To show that these corrections are negligible for the range of 
times s and t where we observe the collapse of the scaling function, we studied the 
two-time correlation C{t,s) on the family of curves tins /tins = k where k is a 
constant. Along these curves, the scaling function gc{i^) is constant and the scaling 
form (5) reduces to s""^ . We then defined an effective exponent as 

lnC{t2,S2) -lnC{ti,Si) 

aes{si, S2) = ^-j — - — ^ (8) 

lns2 — msi 
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Figure 6: Scaling function 8°"= C {t, s) of the spin-spin autocorrelation function versus 
t/s. The trend does not show any real trend to scaling. 




Figure 7: Scaling function s°''=C{t, s) of the spin-spin autocorrelation function versus 
tlns/slnt. A perfect collapse of the curve for the three largest waiting times can 
be observed in the inset. 
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sin t /tins 

Figure 8: Effective exponent Oeff estimated by a fit of the correlation function as 
C{t,s) ~ s~°"'. The time t is chosen such that the scahng variable t\ns/slnt (or 
t/s in the inset), and thus the scaling function gc in Eq. 7 {gc in Eq. 6), is kept 
constant. The legend indicates the two waiting times (si, S2) used for the calculation 
of the effective exponent. 

where Insi/ti Insi = ^2 In 52/^2 In S2 = In practice, to allow for non-integer 
values of ti and t2 the logarithm of the autocorrelation function C{t, s) was in- 
terpolated linearly between the two nearest simulation points, i.e. hiC{t,s) = 
{t - Int t) lnC(Int t, s) + {1 - t + Int t) lnC(Int t + l,s) where Int t is the largest 
integer smaller or equal to t. 



The effective exponent is expected to tend to the exact value Uc = 0.25 in the 
limit of large waiting times s. At finite values of s, corrections to the behavior 
s~"'= should lead to effective exponents that differ from Uc = 0.25. Our data arc 
plotted on Figure 8. One can see that the effective exponent displays as expected 
a plateau around the exact value Oc = 0.25 when sufficiently large waiting times 
s > 160 arc used. The same analysis have been performed along the curve t/s = 
K = Cst. It leads to an effective exponent significantly larger than the exact value 
ttc = 0.25. Moreover the absence of any plateau of the effective exponent excludes 
the possibility of a value of Uc different from 0.25. This analysis confirms that the 
scaling form that is the most compatible with our data is (7) and not (6). 

3.3 Determination of the autocorrelation exponent Ac 

We are now interested in the autocorrelation exponent Ac defined by the assumption 

that the scaling function gc decays asymptotically as gc{x) ^ ^""^^/^^ where x = 
f In s/s Int is the scaling variable. To take into account the possibility of corrections 
to this behavior, we studied the effective exponent {\/ Zc)eS.{^mm,x^a.yL) defined 
as the decay exponent of the autocorrelation fimction C{t,s) with x = t In s/s Int 
in the window x S [a;min; a:max]- The upper bound Xmax is the largest scaling vari- 
able allowed by our data, i.e. Xmax = x{s,t = 10000). The low bound of the 
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Figure 9: Effective exponent \c/zc obtained by a power-law fit C{t,s) 



~ X 



where x is the scahng variable, i.e. x = tins/ slnt (main figure) when taking into 
account topological defects and x = t/s (insert) otherwise. The fit is performed in 
the range x e [xmin, a^max] where Xmax is the maximum value of the scaling variable 
allowed by our data. Xnun is varied and the effective exponent Xc/zc is plotted 
with respect to 1/xmin- Since we are interested in the asymptotic regime where 
corrections to scaling can be neglected, the interesting region corresponds to a;inin 
large. 



interpolation range is varied between x{s,s) = 1 and Xmax- The data are plotted 
with respect to l/x^i^ on Figure 9. For small values of Xmin, almost all numeri- 
cal data, including small values of x, are included in the power-law interpolation. 
As a consequence, the corrections may play an important role and lead to an ef- 
fective exponent possibly quite different from the asymptotic value \c/ Zc- When 
Xiain is increased, the effective exponent is expected to tend to Xc/zc- When Xmin 
approaches Xmax, the number of points that contribute to the fit decreases and the 
effective exponent becomes noisier. The exponent Xc/zc has to be estimated in an 
intermediate regime. 



The data are plotted on Figure 9. For large values of Xmin, the effective ex- 
ponent stays inside a range corresponding to Xc/zc ~ 1.02(2). This means that 
the autocorrelation exponent is Ac ~ 2.04(5). and therefore appears to be satu- 
rated, i.e. Ac = d. The ratio Xc/zc remains compatible with our estimate 1.02(2) 
when considering more complex scaling variables involving logarithmic corrections 
as X = t{\ns)^ /s{\nt}^ or x = tln-^/sln-^. The same procedure with the scaling 
variable x = t/s leads to an incompatible estimate (Ac/2;c)eff — 0.91(2). This dif- 
ference may be due to the fact that the asymptotic regime is seen only for larger 
times if one does not take into account logarithmic corrections. Note that the de- 
cay exponent 0.86 was obtained for the AFIT when no logarithmic correction is 
considered [5]. One can thus assume that taking into account these corrections, an 
exponent closer to 1 would have been obtained. 
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4 Conclusions 



We have given numerical evidences in favor of the existence of topological defects in 
the paramagnetic phase of the FFIM. The dynamical exponent at the critical tem- 
perature Tc = is estimated to be 2; = 2.004(9), in excellent agreement with the 
value z = 2 expected for the AFIT in the scenario involving topological defects [7] 
and fully incompatible with the Monte Carlo study where the signature of these 
topological defects was not seen. We have then studied the decay of the spin-spin 
autocorrelation function C{t,s). Assuming that C{t,s) scales as the homogeneous 
XY fcrromagnct, a good collapse of the scaling function is observed for large waiting 
times. Moreover, the exponent ac is correctly obtained. This would not be the case 
if logarithmic corrections were not taken into account. Interestingly, the autocor- 
relation exponent is compatible with the simple value Ac = d. Note that it is also 
the case for the Ising chain whose critical temperature is = like the FFIM. 

An important question remains unsolved: in the AFIT, the topological defects 
are known to be vortices. The couplings are not homogeneous in the FFIM so it is 
not obvious to us how these topological defects look like for this model. Unfortu- 
nately, they cannot be observed on snapshots of the spin configurations. Note that 
it is also the case for the AFIT. Moore ct al. managed to pin vortices and make it 
observable by averaging over a large number of spin configurations. Without an idea 
of the shape of these topological defects in the FFIM, we have been unable to apply 
this procedure to the FFIM. Further investigation in this direction is necessary. 
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